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A technique is presented for initializing multiple discrete finite element model (FEM) mode sets for cer- 
tain types of flight dynamics formulations that rely on superposition of orthogonal modes for modeling the 
elastic response. Such approaches are commonly used for modeling launch vehicle dynamics, and challenges 
arise due to the rapidly time- varying nature of the rigid-body and elastic characteristics. By way of an energy 
argument, a quadratic inequality constrained least squares (LSQI) algorithm is employed to effect a smooth 
transition from one set of FEM eigenvectors to another with no requirement that the models be of similar 
dimension or that the eigenvectors be correlated in any particular way. The physically unrealistic and con- 
troversial method of eigenvector interpolation is completely avoided, and the discrete solution approximates 
that of the continuously varying system. The real-time computational burden is shown to be negligible due to 
convenient features of the solution method. Simulation results are presented, and applications to staging and 
other discontinuous mass changes are discussed. 


I. Introduction 

In integrated flight dynamics formulations for launch vehicles, the most common approach to modeling elasticity is 
the superposition of a set of linear orthogonal elastic modes with the nonlinear rigid body dynamics. This type of model 
yields many convenient benefits for control design and stability analysis 1,2 . The elastic modes are generated from a 
finite element eigensolution and the corresponding rigid body modes are removed via partitioning. This technique can 
yield a very high fidelity simulation if the nonlinear model elements and the elastic elements are carefully coupled 
together using mathematically rigorous mechanisms. However, the generation of an eigensolution is fundamentally a 
linear, time-invariant (LTI) concept and does not transfer seamlessly to modeling the elasticity of systems that have 
time- varying mass and stiffness, as is common in launch vehicles. 

Alternative methods of modeling require a different formulation and often highly increased computational com- 
plexity. The method of assumed modes 3 fits elasticity to a linear combination of a set of basis functions that are not 
necessarily orthogonal to the rigid body, so the coupling between the elastic response and the rigid body motion has 
to be explicitly applied 4 . Determination of the time- varying terms requires at least partial knowledge of the mass and 
stiffness properties of the original FEM, and increases the computational burden significantly. However, the underly- 
ing structural properties (mass and stiffness) can be interpolated, yielding a quasi-continuous elastic model that closely 
resembles a time-varying physical system. This approach is well-suited to certain types of simulation environments, 
but it is not particularly advantageous for GN&C flight dynamics simulations. 

In a mature launch vehicle design environment, it is typical to have a separate finite element model corresponding 
to a different vehicle mass at every ten to twenty seconds of flight time. At those mass configurations where the sim- 
ulated (table-lookup) rigid body mass properties match the mass properties of the finite element model that generated 
the eigensolution, the elastic motion is indeed orthogonal to rigid body motion. The exact degree of this error can be 
computed given sufficient data from the FEM, but it is known to be negligibly small over the time scales typically con- 
sidered in a launch vehicle dynamics model. When each new model is selected, some mechanism must be employed 
to transition from integrating the states of the eigensolution at time to to the new eigensolution at time t \ . For simple 
vehicle configurations whose elastic dynamics closely resemble that of a beam, it is common to attempt to linearly 
interpolate the mode eigenvectors. Unfortunately, there is no mathematical basis for doing so; the concept of mode 
interpolation is based on the spatial concept of an eigenvector as a representation of the lateral elastic deformation 
of a beam. Mode interpolation requires not only interpolation between similar eigenvectors but also the stiffness and 
damping matrices that represent the elastic generalized coordinates, introducing several problems. 

1 . For continuity, the elastic models must be of the same dimension k at time to and t \ . 
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2. Interpolation between any pair of eigenvectors vq and v\ requires that they are a similar pair; that is, representing 
the same physical motion. This requires a mechanism to correlate modes with physical motion as the models 
evolve forward in time. Automation of this process using MAC (Modal Assurance Criterion) and other tech- 
niques has thus far been elusive, and is highly labor-intensive. This issue is exacerbated by the eigenvectors’ 
arbitrary rotation with respect to the structural coordinate frame. The process does not scale to multiple models 
(such as dispersed model sets), since each model must be preprocessed extensively. For complex models with 
three-dimensional, non-beam-like behavior, this process may not yield a tractable solution at all. 

3. Interpolation of diagonal generalized stiffness and damping matrices is physically intuitive, but the integrated 
model stiffness and damping may no longer be diagonal or symmetric when augmented with generalized forces 
such as aeroelasticity 5 . 

4. Interpolation of several multi-dimensional tables at every integration time step is computationally expensive. 

5. The interpolated model, rigorously speaking, has no correlation to a physical system. 

6. The interpolated model not only non-orthogonal to the rigid body, it is non- orthogonal to itself to an unknown 
degree! 

Even when using modal interpolation with limited success, the problem still remains as to how to transition to a new 
elastic model when the dimension of the model changes, or the model changes in a way that does not lend itself 
to interpolation. Such events include staging or other significant losses of mass that drastically change the modal 
solution. Transfer of the terminal conditions at the last integration time step to the initial conditions of the new model 
is obviously problematic - excepting dimension changes, the new model configuration is represented by an entirely 
new basis, so the old eigenstates are physically meaningless. Continuing the integration often results in a severe, 
physically unrealistic “twang” as the lightly-damped structural modes respond to a large, arbitrary initial condition. 

One obvious solution is to re-initialize the initial conditions of the generalized coordinates to zero. This fails to 
capture dynamically significant events resulting from the release of stored energy in the loaded structure, and can be 
considered underconservative. Other techniques have been proposed to effect a transfer of the initial conditions of 
one model to another, including least squares or pseudoinverse solutions. The problem to be solved is thus: given 
a known set of physical bending displacements and velocities at time to, what set of new generalized coordinates 
optimally approximates the physical coordinates using a new basis (the new set of eigenvectors) at t\ ? The problem 
can be solved quite accurately using a sufficient number of eigenvector components (gridpoints) so that the elastic 
displacements and velocities very closely match on either side of the transition. 

However, as will be shown, this solution creates significant problems due to the lack of continuity of energy - a 
least-squares or weighted least-squares solution is insufficient to ensure a physically meaningful transition. A special 
form of least squares, called regularized least squares or least squares with quadratic inequality constraints (LSQI), 
can be employed instead. 


II. Concept Review 


Consider the linear, second-order structural dynamics model 

Mq + ( Kq = Q 


( 1 ) 


where q eW 1 is a vector of displacements (including small rotations) in physical coordinates, and Q e R n is a vector 
of excitation forces and torques. 

Without loss of generality, assume that the system (1) can be partitioned such that 
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where v are the translational displacements and 6 are the angular displacements, M is the mass matrix, J is the inertia 
matrix, and K\,K 2 are the translational and angular stiffness matrices, respectively. For this second order system (2) 
with M, J , K t real and symmetric, there exists a similarity transformation 
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( 3 ) 


that diagonalizes the system. The transformation S is orthonormal and S 1 * 3 * = S T . Letting x = and 6 = ¥ cr and 
premultiplying by S , we have 
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Typically, O, T 1 are chosen such that the generalized mass matrix is normalized, and 
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A fixed viscous damping is usually assumed, leading to a diagonal damping matrix D. The modal elastic model is 
written as 

jaf] + D?) + 51^7 = S t Q. (7) 

The columns of the transformation S are the eigenvectors of the linear transformation that diagonalizes the system; 
the eigenvalues appear on the diagonal of 51 and are equal to the square of the natural frequencies of the elastic modes. 

The kinetic and potential energies of this class of dynamic system are invariant under coordinate transformations. 
It follows, therefore, that the kinetic energy of the system 

T = q T Mq - rfq. ( 8 ) 

Likewise, the potential, or elastic strain energy, is given by 

V = q Tr Kq = q T Siq. (9) 


III. Model Transition Using Least Squares 

Consider two finite element models (1) in modal form (7) at times to and t\ . Given any complete eigenstate 770 , fjo 
at time to, the physical elastic displacement and velocity of the model are given by 

qo = S oqo (10) 

<7o = S o?7o 

Recall that at t\, q\ = S \q\ and likewise for the velocities. We seek a solution for the new eigenstate fji at 
time t\ that closely approximates the physical displacement qo such that 1/71 - qo\\ 2 , the norm error of the approxi- 
mation, is minimized. This is a standard estimation problem and a solution is readily found using the Moore-Penrose 
pseudoinverse given by 


7/i = (ii) 

m = 

where 

5j = (5[5i)"5[ (12) 

which is guaranteed to exist since S \ is full rank. " 

Although the pseudoinverse solution ensures the minimum square error of the physical approximation, it does so 
with arbitrary 77 , 77 . If it were the case that the model (7) contained the full set of modes corresponding to all of the 
modeled mass in the original FEM, the matched displacement and velocity would provide approximately the same 
strain energy and kinetic energy, respectively, at time t \ . However, the FEM (7) is often significantly truncated, so the 
physical solution is an approximation, S t are not square, and energy continuity is not guaranteed and is in fact grossly 
in error. 

'The matrix 1 can be very ill-conditioned; direct use of (12) may yield errors. More numerically robust methods for computation of the 

pseudoinverse are common; for example, see MATLAB’s “\” and pinv() functions, which employ QR decomposition 6 . 
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A. Solution Via Regularized Least Squares 

A least squares problem is said to be regularized when the model parameter estimate x eR n that achieves 

min\\Ax-b\\ 2 (13) 

is subject to lie within a hyperellipsoid in R”, which can be defined as the norm inequality 

\\B(x - xo )||2 < a. (14) 

The regularized problem is useful in certain estimation applications, for example, where a nonlinear process is 
linearized at xo and a model parameter estimate is refined via iteration. At each step, the new least-squares estimate 
for the parameters is constrained to lie within a trust region around the current linearization. 

There exist algorithms 6,7 that readily solve the regularized least squares problem (13,14), relying on a generalized 
singular value decomposition (GSVD) which diagonalizes the problem into a set of scalar objectives and constraints. 
The method requires the determination of a Lagrange multiplier via Newton-Rhapson iteration. 

Note the equivalence of the quadratic inequality constraint 

(x - x 0 ) T Q{x- v 0 ) < a 2 (15) 

where Q = B T B, a symmetric, positive definite matrix. Given any symmetric, positive-definite (and therefore diago- 
nalizable) Q , we can always find a B and write the inequality in the form of (14). The quadratic form in (15) can be 
used to express the potential and kinetic energies (8) and (9). We reformulate the two least squares problems as 

min \\S ir/i - qo\\ 2 subjectto |^ 57 7 i| 2 ^ V^o (16) 

which determines the generalized displacements such that the new potential energy V\ < Vo- For the velocities, we 
find 

min IISi^i -<7oll 2 subject to ||t 7 i|| 2 < (17) 

such that the new kinetic energy T\ < T 0 . Note that as a consequence of the aforementioned diagonalization, V 0 = 
tIqRtio and To = rjodo- Furthermore, the square root of the matrix required in (16) is trivial since 51 is diagonal. 

B. Numerical Implementation 

Numerical implementation of the LSQI constrained least squares problem (16,17) is straightforward using MATLAB’s 
built-in SVD and matrix routines along with the open source (BSD license) package regtools* , which can be readily 
adapted to the present application. Of the three distinct constraint cases occurring in the LSQI problem, we find that 
the numerical solution invariably lies on the boundary of the feasible set 6 ; that is, the best-fit solution for the physical 
displacements and velocities is the maximum-energy solution. 

In practical implementation, it may not be necessary to include the entire set of eigenvectors Si but only the 
translational degrees of freedom associated with a suitable number of gridpoints distributed throughout the structure 
such that global elastic behavior is well-represented. The number of required gridpoints should equal or exceed 
the number of generalized degrees of freedom, but varies depending on the degree of model truncation. Adequate 
fitment with only low-frequency global bending modes can be accomplished with only a few dozen gridpoints. The 
elastic displacement and velocity of the gridpoints used for initialization need only be computed when executing the 
algorithm. There is no requirement that the models at to and t\ be of the same dimension. 

More importantly, the numerically intensive generalized singular value decomposition (GSVD) can be computed 
offline based on the modal data; the result is not dependent on the current states. Therefore, the initialization technique 
is well- suited to simulation implementation since the GSVD matrices associated with each model can be stored in a 
lookup table. 


IV. Simulation 

The present initialization technique has been applied to a high fidelity, linear time- varying (LTV) planar perturba- 
tion simulation of the Ares I launch vehicle derived from the FRACTAL 9 modeling toolchain. In order to adequately 
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model elasticity effects on control system performance, the model can be initialized to propagate approximately 16-30 
modes; global bending behavior can be captured sufficiently with 8-10 modes. 

In the present case, 16 lateral bending modes have been down-selected from the elastic model using a modal gain 
sensitivity metric 10 . The initialization routine uses 53 centerline gridpoints along the vehicle core load path. No 
frequency interpolation is used. Since a consistent model dimension is used, the number of modes at each transition 
is the same. The transition algorithm is implemented in MATLAB using a modified version of regtools. Under these 
conditions, the transition algorithm typically converges in about 2 ms on a 3.0 GHz Intel E8400 PC platform and 
therefore has no impact to overall simulation performance. 

The test cases considered are a nominal ascent simulation with multiple flex model transitions; the first occurs at 
t = 5 seconds and the second occurs at t - 15 seconds. The second transition is analyzed to demonstrate the efficacy of 
various transition methods. The blended rate output sensed by the control system (consisting of sensor measurements 
from both forward and aft rate gyro assemblies) yields a qualitative measure of the amount of flexibility affecting 
the autopilot. In addition, the elastic potential V(t ) is computed. In the following plots, the elastic potential has been 
forward-reverse processed with a 4th order Butterworth low-pass filter (f c = 2 Hz). The simulation integration interval 
is 0.01 seconds using an RK4 algorithm. 

Three test cases are considered. First, the transition to a new set of eigenvectors is allowed to occur without a 
change in the generalized coordinates (no initialization). The second test case initializes the new generalized coordi- 
nates to minimize the weighted square error of the physical displacement and velocity, as in (11). Finally, the third 
test case uses the FSQI energy-based method (16,17). The method of resetting the generalized coordinates to zero is 
not shown; the results are obvious. 

In Figures (1) and (2), the transition is allowed to occur without initialization. No significant increase in bending 
response at the sensor can be detected, and energy continuity is approximately maintained. Since the generalized 
coordinates do not change, the kinetic energies are equal and the potential energy is nearly the same. However, as 
can be seen in Fig. 2, the physical configuration after the transition is in no way related to the physical configuration 
before. 



Figure 1. Elastic angular response and potential energy, no state initialization 
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Flex init at T=14.990 




Figure 2. Mode shape and velocity continuity, no state initialization 


A standard pseudoinverse initialization using the minimum square error of the displacement and velocity is shown 
in Figures (3) and (4). Due to the ability of the new basis to accurately model the preceding physical configuration, 
the model fit in physical coordinates is very good. However, the problem with the unconstrained solution becomes 
apparent - the bending potential discontinuously increases by several orders of magnitude, creating a large elastic 
transient that rapidly drives the numerical simulation to instability. 



Figure 3. Elastic angular response and potential energy, pseudoinverse initialization 
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Flex init at T=14.990 




Figure 4. Mode shape and velocity continuity, pseudoinverse initialization 


Finally, the proposed LSQI technique is illustrated in Figures (5) and (6). Kinetic and potential energy continuity 
is maintained across the transition, and a constrained solution exists that very accurately approximates the physical 
configuration of the system before the change of basis vectors. No unrealistic transient is introduced into the system. 
The simulated elastic behavior through the transition closely matches the expected physical elastic response. 



Figure 5. Elastic angular response and potential energy, LSQI initialization 
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Flex init at T=14.990 




node number 

Figure 6. Mode shape and velocity continuity, LSQI initialization 


V. Conclusions and Forward Work 

The proposed LSQI initialization technique has been shown to provide physically accurate simulation continuity 
for differing elastic models without the requirement that the models be of the same dimension or have similar sets 
of eigenvectors. This algorithm may yield sufficient reductions in computational and manual preprocessing overhead 
to allow the use of multiple sets of dispersed modal solutions, since manual or automated inspection of mode shape 
correlation can be eliminated. More importantly, the speed and fidelity of GN&C flight dynamics simulations can be 
increased by eliminating unnecessary table interpolation steps and improving the physical realism of model transitions 
while retaining the simplicity of the orthogonal modes approach. 

The sensitivity to the number of modes and the number of gridpoints used in the transition algorithm remains to be 
shown. Convergence of the embedded iterative root finding step is not guaranteed, and it may be possible to improve 
convergence properties by including the rotational degrees of freedom in the transition algorithm. Convergence might 
also be improved in a second pass by re-initializing the LSQI solver using a random guess for the new eigenstate that 
lies on the boundary of the constraint ellipsoid. 

One particularly intriguing applications area is enforcement of continuity during discontinuous changes in the 
vehicle mass or configuration, such as a panel jettison or staging event. This requires knowledge of the fraction of 
kinetic and potential energy that is retained after the transition so that the remaining structure to be propagated forward 
is properly initialized. 
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SPACECn, 


Overview 



♦ Many 6+ DoF dynamics formulations are based on superposition of linear 
orthogonal elastic modes with nonlinear rigid body dynamics 

♦ A FEM eigensolution (set of modes) is fundamentally an LTI concept 

♦ The vehicle dynamics are not LTI; system has time-varying mass/stiffness 

♦ Superposition of the orthogonal modes is still the most attractive modeling 
approach 

• Does not require knowledge of mass integrals; coupling between rigid and elastic DoF 
found in assumed modes techniques 

♦ In a mature modeling environment, FEMs are expected at closely spaced 
intervals 


• Models are only orthogonal to rigid body at these distinct mass property configurations 

• Error is small over these time scales using fixed modal solution 

• How do we transition from one model to the next? 
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Mode Interpolation 

(OR, WHY IT’S SUCH A BAD IDEA) 


♦ Mode interpolation has been used in the community for some time 

• Mode interpolation is physically intuitive for simple beam systems 

• For general structural dynamics solutions, mode interpolation does not yield a physically 
meaningful answer 

1 . For continuity in integration, model dimension must not change 

2. Interpolation between eigenvectors requires that they be a similar pair 

• Introduces significant preprocessing overhead (MAC, human analyst) 

• Eigenvectors can be arbitrarily rotated w.r.t. structural frame (“sign flipping ’) 

• Does not scale to sets of dispersed models 

3. Inclusion of aeroelastic & force follower effects makes stiffness/damping non-diagonal 

4. Multi-dimensional table interpolation is computationally expensive 

5. Interpolated model, strictly speaking, is not correlated to any physical system 

6. The interpolated model is non-orthogonal to itself 
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Further Complications 


♦ Elastic model transitions at discontinuous events 

• Discontinuous events require propagation of a new elastic model with possibly no dynamic 
similarity to previous model 

• Consider loss of mass (staging), panel separation, LAS jettison, etc. 

• Dimension of the model may change 

♦ Existing solutions are either unrealistic or underconservative 

1 . Transfer initial conditions (or subset of them) to the new model 

No physical correlation to the real system; represented by an arbitrary new basis 
Analogous to expressing a physical 3-vector in the wrong coordinate frame 

2. Transfer initial conditions and artificially increase modal damping to absorb twang 

Not physically realistic; fixes a modeling problem 

3. Set initial conditions to zero 

Underconservative; ignores or resets loaded structure 
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Example 


♦ Example (Ares I FRACTAL LTV, 16+ modes) 

• Direct transfer of initial conditions to new model at T=1 5 sec 

• No significant increase in energy (kinetic energies equal; small AV due to Aw) 

• No physical correlation between old and new system representation 



Elastic Potential Energy 



Flex init at T=1 4.990 
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Concept Review 


♦ Linear, 2 nd -order structural dynamics model 


Mq + Kq = Q q e R n Q G M n 


• Assume we can partition this model as 
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Concept Review (II) 


♦ We write the structural model in a modal form 

f-irf + D?7 + J li] = S r Q 

• The kinetic energy is invariant w.r.t coordinate transformations - 

T = q T Mq = V T q 

• As is the potential (elastic strain energy) 

V = q T ICq = ifsb) 

• Note that in modal space, the energy can be expressed as a quadratic form of a vector and a 
diagonal (or identity) matrix 



Model Transition Using Least Squares 


♦ How do we effect a smooth model transition? 

• Given a known elastic displacement using model “A” (set of basis functions [eigenvectors]), how do 
we find a new set of generalized coordinates that make “B” displacement match “A”? 

♦ Consider two finite element models in modal form at and ^1 

• Elastic displacement and velocity given by 

g 0 = Sqtjq 
Qo — SoVo 



• Minimize |Si?7i — ®)|| 2 : 


m 

= si go 

m 

= S'l 9o 


• This is a least-squares model fit; a Moore-Penrose pseudoinverse: 



(sfs^sl 


• Attractive and intuitive solution; but it doesn’t work! No energy continuity. 

• Physical displacement and velocities match very well, but only at time of transition. 

• This is partially due to the truncated model; modal mass in each model is not the same. 
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Model Transition Using Least Squares 


♦ Example (Ares I LTV, 16+ modes) 


Flex init at T=1 4.990 
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Model Transition Using Least Squares 


♦ Example (Ares I FRACTAL LTV, 16+ modes) 

• Unstable numerical response after transition due to large increase in elastic potential 


Rate Channel Elastic Response 
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Model Transition Using LSQI 


♦ Consider the previous problem as a regularized least squares problem 

• In estimation theory, we regularize a least squares problem to confine the new estimate to be close 
to the old estimate, for example in successive estimation of a nonlinear process 

• The least squares procedure minimizes 


• Subject to 


min | Ax — b\ 0 
||S (x — ) 1 1 2 


• B (its square, actually) defines a hyperellipsoid in 

• If B=l, x must lie in a unit ball around 


• There are existing algorithms that solve this nonlinear problem efficiently. 


♦ Note the equivalence of the quadratic inequality constraint 

(z - zof Q i x ~ x o) < a 2 Q = B T B 

and our forms for the kinetic and potential energies, 

T — q T Mq = r} T i} 

V = q T )Cq = n T & ; 
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Model Transition Using LSQI (Cont’d) 


♦ We now solve the quadratic inequality constrained least squares (LSQI) problem 


Match displacement subject to potential energy constraint 


min | jSi 771 


<70 II2 subject to 




<VVo ■=> U < v 0 


Match velocity subject to kinetic energy constraint 


min | Si 7/1 — qq 1 1 2 subj ect to ||/)i ll 2 < cz> < To 

♦ Since modal form is diagonal, quantities are easy to compute: 

V 0 = rjo&r]o T 0 

♦ Numerical implementation of LSQI is available in an open-source MATLAB package, 
regtools, modified slightly for this application 

• Only a small number of centerline gridpoints need be computed, only at the model transition 

• The GSVD (generalized singular value decomposition) calls are not a function of the current states: 
all complex matrix math can be computed offline 

• Numerical routine requires a root-finding step to determine a Lagrange multiplier 

• Considerably less computationally complex than interpolation 
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Simulation 



♦ Flex initialization routine implemented in Ares I FRACTAL LTV, 16+ modes 

• 53 centerline (load path) gridpoints used for displacement/velocity matching 

• 200 Hz continuous-time integration rate on analyst PC (3.0 GHz Intel E8400) 

• Full aeroelastics, thrust, and slosh coupling models 

• Typical execution time (including GSVD) < 2 ms! 

• Eliminates twang in a physically realistic way 
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Conclusions / Forward Work 



♦ LSQI algorithm can provide physically accurate simulation continuity without modal 
interpolation 

• Only losses of fidelity are (1) modal truncation (2) loss of orthogonality w.r.t. rigid body 

• Ability to match displacement and velocity of gridpoints at separation events to simulate pre-load 
induced twang effects 

- Could also be used to transition from fixed-free to free-free stack model at liftoff 

- If mass loss occurs, knowledge of remaining mass/stiffness distribution in each “piece” can be 
used to determine new energy constraints for separate free-free bodies 


♦ Potential to significantly reduce preprocessing overhead and computational burden 

• Table lookup of interpolated modal data is very expensive 

• Manual preprocessing of large, complex three-dimensional structural models is nearly intractable 

• Allows automated use of dispersed (realistic) modal models in GN&C simulations 


♦ Forward work 

• There is an embedded root-finding problem whose convergence properties must be explored 

• Analysis of sensitivity to number of modes and accuracy of model fitment must be explored 

• Inclusion of rotational degrees of freedom may add fidelity 

• Verification of accuracy of discretized solution against continuously varying analytical solutions 
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